Ab Initio Many-body Study of Cobalt Adatoms Adsorbed on Graphene 
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Many recent calculations have been performed to study a Co atom adsorbed on graphene, with significantly 
varying results on the nature of the bonding. We use auxiliary-field quantum Monte Carlo (AFQMC) and a 
size-correction embedding scheme to accurately calculate the binding energy of Co on graphene. We find that 
as a function of the distance h between the Co atom and the six-fold hollow site, there are three distinct ground 
states corresponding to three electronic configurations of the Co atom. Two of these states provide binding 
and exhibit a double-well feature with nearly equal binding energy of 0.4 eV at h — 1.51 and h = 1.65 A, 
corresponding to low-spin 2 Co (3rf 9 4s°) and high-spin 4 Co (3d 8 4s 1 ), respectively. 

PACS numbers: 61.48.Gh 73.22.Pr 73.20.Hb 31.15.A- 



Since its discovery, graphene has been the subject of in- 
tense efforts to adapt it for a variety of promising applica- 
tions due to its unique and exceptional intrinsic properties J£ 
One potential application is for use in spintronic devices." 
However, external methods are required to induce magnetism 
on graphene, since pristine graphene is nonmagnetic. One 
proposal is to adsorb transition metal atoms to provide lo- 
calized magnetic moments in graphene. Single Co atoms on 
graphene have been extensively studied recentlyr^ and pos- 
sible Kondo effects have been consideredi22i2i The study of 
Co/graphene is thus of great interest both from a fundamental 
and applied perspective. 

Theoretical treatments of Co/graphene systems have largely 
been done at the density functional theory (DFT) level with 
local or semi-local functionals, or with an empirical Hubbard 
on-site repulsion U (DFT+U)&~~ However, the applicabil- 
ity of methods based on independent-electron approximations 
in such systems is unclear, since electron correlation effects 
can be significant. Indeed, widely varying results have been 
reported for the nature of the magnetic state and binding of 
Co as a function of adsorption height. DFT calculations with 
the generalized gradient approximation (GGA) 2 ^ predict^— 
an equilibrium height of h eq ~ 1.5 A above the six-fold hol- 
low site, with a low-spin Co atom configuration (S = 1/2). A 
different functional, the hybrid B3LYP*22 predicts^ an equi- 
librium height of h eq ~ 1.9 A at the hollow site, with a high- 
spin configuration (S = 3/2). Results from the GGA+U ap- 
proach have shown sensitivity to the choice of the parame- 
ter U which leads to different spin configuration, equilibrium 
height, and equilibrium site for different values of U i 10 ' 14 ' 18 A 
recent quantum chemistry calculation using the complete ac- 
tive space self-consistent field method gives a state from the 
van der Waals (vdW) interaction (high-spin 3d 7 4s 2 state) as 
the global minimum, with /i eq ^3.1 A^ These contrasting 
results strongly indicate the need for a more accurate ab initio 
treatment of electron correlations in Co/graphene. 

In this paper, we use the auxiliary-field quantum Monte 
Carlo (AFQMC) metho d 24 ' 25 to investigate the binding energy 
and electronic properties of Co/graphene. We focus on the 
hollow site which is the most favorable adsorption site accord- 
ing to most DFT calculations. Contrary to prior calculations, 
we find that as the Co atom approaches the graphene sheet, 
it experiences two magnetic transitions which lead to three 



distinct ground-state electronic configurations. One of these 
configurations corresponds to the vdW interaction. The other 
two configurations arise from a strong orbital hybridization 
and provide binding with a double-well feature. 

Since strong electron-electron interactions are expected to 
be spatially localized in the immediate vicinity of the Co atom, 
we use a size-correction embedding scheme (ONIOM 2 ^) to 
accelerate convergence and reach large system sizes in the 
many-body calculations. In this approach, the "near" region 
in the vicinity of the Co atom is modeled by a relatively small 
number of atoms, using a highly accurate many-body method 
like AFQMC, while size corrections are treated using a lower 
level of theory like DFT. For the near region, we chose the 
Co atom and its six nearest neighbor substrate C atoms, with 
their dangling bonds terminated by H atoms, resulting in a 
Co/CgH6 benzene-like system (see the inset in Fig. [T). The 
size-corrected binding energy of the Co/graphene system is 
then given by 



ONIOM 



= E 



Co/C 6 H 6 
b, AFQMC 



/ j-,Co/graphene 
l^b, DFT 



b, DFT )> 



which we will calculate as a function of h, the perpendicu- 
lar distance between Co atom and the substrate, for each spin 
multiplicity of the Co atom. For each substrate, E\, is defined 
as E h = £ Co/substrate - E Co - E suhstrRte . The Co/C 6 H 6 C- 
C bond length was fixed to that of graphene, 1.42 A, which 
is only slightly larger than the experimental benzene value 
of 1.40 A, while the distance to the H "link atom," the C-H 
bond length, was set to 1.09 A, which is the predicted geom- 
etry by GGA for the corresponding C-C bond length. Previ- 
ous studies have shown little sensitivity to the link-atom bond 
distance. 27 Our AFQMC calculations were all done for fixed 
substrate geometries. We will consider the effect of substrate 
geometry relaxation with the assistance of DFT calculations, 
as discussed below. 

The AFQMC method^l^ 5 . evaluates the ground state proper- 
ties of a many-body Hamiltonian stochastically, using random 
walks with Slater determinants expressed in a chosen one- 
particle basis. Although AFQMC is an exact method in princi- 
ple, the fermion sign problem causes an exponential growth of 
the Monte Carlo variance. The problem is controlled using a 
constraint on the overall phase of the Slater determinants dur- 
ing the random walks, the phaseless approximation, 25 that re- 
lies on a trial wave function. In extensive benchmarks in both 
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FIG. 1. (Color online) Binding energy of Co on CqHq as a func- 
tion of Co adsorption height h at the six-fold site for different meth- 
ods. For AFQMC, left, middle, and right curves correspond to nom- 
inal 3d 9 4s°, 3d 8 4s 1 , and 3d 7 4s 2 Co configurations, respectively. 
AFQMC results include Trotter time step extrapolation. For DFT re- 
sults, the left and right curves correspond to 3rf 9 4s° and 3ci 8 4s 1 Co 
configurations, respectively. The shaded area on the AFQMC Morse 
fits reflects one standard deviation statistical errors. 



strongly correlated lattice models and molecular and crys- 
talline systems, the method has shown excellent agreement 
with exact and/or experimental results i 24,25 ' 28- — This is con- 
sistent with expectations from analysis of the origin of the sign 
problem and the nature of the constraint, 24 ' 25 In most calcula- 
tions to date on realistic systems (molecules and solids), trial 
wave functions of a single Slater determinant from Hartree- 
Fock or DFT have been used and have been shown to give 
results whose accuracy is comparable to the best many-body 
methods, for example coupled-cluster CCSD(T) in molecules. 
In this paper, we use the phaseless AFQMC method working 
with standard Gaussian single-particle basis sets* 2 ^ and a re- 
cent implementation of the frozen-core approximation to treat 
the inner core electrons.— 

We first report AFQMC results for Co on a CqHq substrate. 
In themselves, these results provide a direct and systematic 
benchmark of other computational methods. The binding en- 
ergy curves of Co/C 6 H 6 from AFQMC and DFT (GGA and 
B3LYP) are shown in Fig.Q] AFQMC results show that, with 
decreasing h, the ground-state electronic configuration of the 
Co atom undergoes two transitions resulting in three differ- 
ent configurations: high-spin id 7 As 2 , high-spin 3rf 8 4s 1 , and 
low-spin 3d 9 4s states, respectively. Only two DFT ground- 
state configurations are found, a high-spin 3(i 8 4s 1 for large h 
and a low-spin 3d 9 4s° for small h. This is because both DFT 
functionals incorrectly predict 3c? 8 4s 1 to be the ground state 
configuration for the free Co atom. Both GGA and B3LYP 
predict low- and high-spin relative minima in the vicinity of 
the AFQMC predictions, but GGA severely overestimates the 
well-depths, which are underestimated by B3LYP. 

The AFQMC calculations were done with our recently 
implemented frozen-core approximation, 33 thus avoiding the 
need for pseudopotentials; only the most tightly bound core 
states were frozen: Co(ls, 2s, 2p) and C(ls). The potential 
energy curves (PECs) are obtained by AFQMC calculations 



TABLE I. Calculated binding energies Eb and adsorption heights h 
of Co on CeHe at the three local minima shown in Fig. Q] (distances 
in A and energies in eV). Tabulated Eb at the low-spin h = 1.47 A 
and high-spin h = 1.65 A minima are CBS extrapolations. Eb at the 
van der Walls (vdW) h = 3.4 A minimum is essentially converged 
at the QZ level. 
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with fixed S z , in which the numbers of electrons with 'f- and 
^-spins are preset. The Co spin configuration in the differ- 
ent PECs is identified by that of the trial wave function 
Thus these are nominal states and do not imply literal spin 
configuration of Co in the many-body ground state. Typical 
AFQMC runs used ~ 5000 walkers and a Trotter time step 
At = 0.01 Ha -1 , and final results were extrapolated to the 
At ->• limit. All AFQMC calculations for Co/C 6 H 6 used 
single-determinant unrestricted Hartree-Fock (UHF) \I>t- Pre- 
vious experience indicates that such AFQMC calculations are 
very accurater^ 3 ^ 4 - including for systems containing transi- 
tion metal atoms. 28 Future study using multi-determinant "fx 
is warranted, however, given the new territory being explored 
here with QMC. DFT and HF calculations that use Gaussian 
basis sets are performed using NWCHEM, 35 - 

Care was taken to remove finite basis set error in the many- 
body results. The following basis sets were used in AFQMC 
calculations for most h. The Co atom used the correlation- 
consistent core-valence cc-pwCVTZ basis set, where "core" 
refers to the Co 3s, 3p semicore states. For C and H atoms, 
valence-only cc-pVTZ and cc-pVDZ were used, respectively. 
For several geometries near the minima, the Co(cc-pwCVQZ) 
basis set was used to obtain extrapolation to the complete ba- 
sis set (CBS) limit. Not surprisingly, while the DFT results 
are converged by the Co(cc-pVTZ) level, AFQMC is not yet 
fully converged even at the Co(cc-pwCVQZ) level. To esti- 
mate the effect of the CBS extrapolation, we used the proce- 
dure in Ref. [3^: an exponential form for the HF contribution 
to the total energy and an inverse-third-power form for the 
correlation energy. Extrapolation to the CBS limit lowers the 
binding energy near the minima by 0.13eV from the TZ re- 
sult and 0.03 eV from that of the QZ basis. Trotter time step 
extrapolations were obtained from results for a smaller basis 
set [cc-pVTZ for Co and cc-pVDZ for C and H] and applied 
to the E\, results for the larger basis sets. Final, fully extrapo- 
lated AFQMC results at the three minima are tabulated in Ta- 
bleUl The global minimum in AFQMC is the low-spin 3d 9 4s° 
state with binding energy Eb = — 1.07(6) eV, as seen in the 
Table. The high-spin minimum has only a slightly smaller 
E\, = — 0.92(5) eV. In the vdW region, the system is barely 
bound with E h = -0.10(3) eV. 

Results are then obtained for Co/graphene using the 



3 




AFQMC-relaxed / cc-pwCVTZ 
AFQMC / cc-pwCVTZ 
AFQMC / cc-pwCVQZ 



1.2 



1.6 



2.0 



2.4 2.1 

h (A) 



3.2 



3.6 



4.0 



FIG. 2. (Color online) Binding energy of Co atom on graphene as a 
function of h. Left, middle, and right curves correspond to 3d 9 4s°, 
3d 8 4s 1 , and 3d 7 4s 2 Co configurations, respectively. Shaded areas 
are one-a statistical error bars. The left inset shows the structure 
of Co on graphene in 5 x 5 supercell. The right inset shows the 
binding energy after CBS extrapolation and substrate relaxation (see 
text). The shaded areas in the right inset include both statistical and 
systematic errors. 




FIG. 3. (Color online) ONIOM size corrections and the binding en- 
ergies of Co/CeHe and C0/C24H12 systems in the 3d 8 4s 1 state. The 
two ONIOM curves are basically identical and show insensitivity to 
the choice of DFT flavors. The corrections are applied to the 3d 8 4s 1 
AFQMC/cc-pwCVTZ binding energy curve in Fig. [2 Similar inde- 
pendence on DFT functional is found for the other spin states. 



ONIOM embedding scheme. The finite-size correction [sec- 
ond term in Eq. ([TJ] is applied to the AFQMC curve in 
Fig- LT] The results are shown in Fig. [2] To obtain E^°l^ hene , 
we used DFT-GGA as implemented in the PWSCF code of 
the QUANTUM ESPRESSO package,. 37 with periodic boundary 
conditions and ultrasoft pseudopotentials.- 8 A 5 x 5 in-plane 
supercell was used, which contains 50 C atoms and a Co 
atom; the in-plane lattice parameter was 12.3 A, while the 
periodic repeat distance perpendicular to the graphene plane 
was set to 15 A. A planewave basis kinetic energy cutoff of 
E cu t = 45 Ry and a charge density cutoff 360 Ry were used 
for all geometries. Brillouin-zone sampling used a T-centered 
4x4x1 fc-point grid and a Gaussian smearing width of 0.04 
eV. The ONIOM E^ correction was obtained from similar 
PWSCF calculations for the clean 5x5 graphene supercell; the 
energy of an isolated Co atom was obtained using a large su- 
percell with single fc-point sampling. Approximate relativis- 
tic corrections are included in our results via ONIOM as the 
GGA calculations are scalar-relativistic, although the correc- 
tion is not perfect due to the absence of the vdW curve in 
GGA. The lines in Fig. [2] are Morse fits to the cc-pwCVTZ 
AFQMC results. 

It is reassuring to note that the size correction in Eq. ((TJ 
is essentially independent of the choice of DFT exchange- 
correlation functional. This is illustrated, for GGA and 
B3LYP, in Fig. [3] using a coronene-like C24H12 substrate 
which is comprised of six joined Cq rings with outer H ter- 
minations. (B3LYP calculations for the 5x5 supercell were 
time consuming and difficult to converge.) As Fig. [3] illus- 
trates, while GGA and B3LYP show large differences between 
their E\, curves, the size correction in Eq. (fTJ is essentially in- 
dependent of which is used. 

We examined substrate relaxation effect by comparing the 
relaxed and unrelaxed 5x5 PWSCF supercell results and in- 
cluding it as an additional ONIOM "layer." For this purpose, 



the six C atoms nearest Co in the relaxed substrates were al- 
lowed to relax only in the in-plane direction. The value of h 
was defined in relation to these atoms; the remainder of the C 
atoms were completely relaxed in C21, symmetry. Relaxation 
was considered complete when the force on all atoms, except 
the restricted atoms, was less than 0.02 eV/A. Near the double 
well minima, fully relaxing all the atoms had little additional 
effect near the low-spin (high-spin) minimum: in- and out-of- 
plane distortions are < 0.015 (0.011) A and < 0.01 (0.002) 
A, respectively. Substrate relaxation lowers the binding en- 
ergy by about 0.05 eV near the minima. 

The inset of Fig. [2] shows the binding energy curves near 
the double well feature, after CBS extrapolation and sub- 
strate relaxation effects have been included. The two wells 
in the Co/graphene PEC's have comparable binding ener- 
gies of — 0.4eV. The vdW region shows no binding within 
AFQMC statistical resolution. STM experiments could, in 
principle, detect the spin-state of Co atoms on graphene.— 
Recently, controllable ionization and screening of Co atoms 
on graphene via scanning tunneling microscopy (STM) have 
been observed^ Kondo screening is generally considered for 
effectively S = 1/2 impurity systems. Higher spin states can 
be observed, however, in the presence of magnetic anisotropy, 
if it results in a low-lying degenerate doublet ground state, as 
was observed for individual Co atoms adsorbed on Cu(100) 
crystals that are covered by a monolayer of copper nitride 
(Cu 2 N).— Mattos^i reported STM observations of Kondo sig- 
natures for Co/graphene. Brar et al.^ however, measured a 
Kondo-like dip feature (with a 5 meV half- width in dl/dV) 
for Co on back-gated graphene/Si02, which they instead at- 
tributed to vibrational inelastic tunneling. To model this 
they performed DFT supercell calculations for free standing 
hollow-site Co/(4 x 4)-graphene and found in-plane vibra- 
tional modes of 12 and 27 meV, and out-of -plane modes of 
17, 40 and 53 meVri the lowest of which are roughly com- 
mensurate with the observed 5 meV width. Within the sta- 
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tistical resolution of the AFQMC double wells in the inset 
of Fig. [2] both low- and high-spin minima have the same 
curvature, corresponding to an out-of-plane frequency range 
16-58 meV, qualitatively similar to the DFT frequencies. 
At liquid He temperatures where the STM experiments are 
performed, tunneling between the minima in Fig.[2]can be ne- 
glected, based on a barrier height of 0.04 eV. Experimental de- 
terminations are further complicated, however, by indications 
that the charge state of single Co atoms on graphene switches 
in proximity to the STM tip.— 1 To the best of our knowledge, 
current experiments for Co/graphene have not yet determined 
the spin state of individual Co atoms adsorbed on graphene. 
Our results are consistent with roughly equal occurrence of 
Co S = 1/2 and 5 = 3/2 atoms populating the two minima, 
respectively. 

In summary, we have presented an ab initio many-body 
study of Co on graphene to address the effect of elec- 
tron correlations. We use the AFQMC method with single- 
determinant trial wave functions to calculate the binding en- 
ergy curve of CoIC§Y{q. The Co/graphene binding energy was 
calculated using an ONIOM size-correction procedure. The 
size-correction method shows insensitivity to the choice of 
DFT flavors which suggests that CoICqH§ cluster captures 
most of the correlation effect. The resulting binding energy 



curve of Co on graphene exhibits binding with a double-well 
structure. Both minima show nearly equal binding energy of 
—0.4 eV. The inner well corresponds to a low-spin 5=1/2 
state with a 3d 9 4s° electronic configuration for Co atom, 
while the outer well is characterized by a high-spin (5 = 3/2) 
3d 8 4s 1 state. Our results show that the Co/graphene system 
requires an accurate and careful treatment of many-body cor- 
relation effects. Better resolution of the energetics and the 
characteristics of the ground states will require further work, 
but the results suggest a plausible framework which is consis- 
tent with recent experimental observations. We hope this re- 
sult will encourage further theoretical and experimental stud- 
ies of the spin states and Kondo effect in Co on graphene. 
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